
/* This program builds up the matrix of forwards based on two sets of initial conditions: (1) an initial forward curve and (2) a time series of the longest-horizon forward rate.
// Matrices are indexed with column major indicies (go down column, then to next column, etc.)

// When calling from matlab, the syntax is:

// [F] = buildForwards(dR, F0, FN)

// F    is the TxN matrix of forwards
// dR   is the TxN matrix of difference returns
// F0   is the NX1 vector of initial forward rates
// FN   is the TX1 vector of the longest-horizon forward rate*/

/********************************************************************************************/

#include <math.h>
#include "mex.h"

/* Define Input Arguments*/

#define Fout_   plhs[0]          /*TxN matrix*/
#define dRin_   prhs[0]          /*TxN matrix (difference returns) [first row = -999]*/
#define F0in_   prhs[1]          /*Nx1 vector (initial forward curve)*/
#define FNin_   prhs[2]          /*Tx1 vector (longest-horizon forward rate)*/

void mexFunction ( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] )

{
    double *dR, *F0, *FN;
    double *F_ptr;
    int T, N;
    int t, n;
    
    N = mxGetN(dRin_);       /*get value of T*/
    T = mxGetM(dRin_);
    
    /*assign pointers to inputs*/
    dR = mxGetPr(dRin_);
    F0 = mxGetPr(F0in_);
    FN = mxGetPr(FNin_);
    
    /*set up output matrix*/
    Fout_ = mxCreateDoubleMatrix(T,N,mxREAL);
    F_ptr = mxGetPr(Fout_);
    
    /*build initial conditions*/
    for (t=0; t<T; ++t) {
        F_ptr[(N-1)*T+t] = FN[t];
    }
    
    for (n=0; n<N-1; ++n) {
        F_ptr[n*T] = F0[n];
    }
    
    /*build forwards*/
    for (n=N-2; n>-1; --n) {
        for (t=T-1; t>0; --t) {
            F_ptr[n*T+t] = F_ptr[(n+1)*T+t-1] - dR[(n+1)*T+t];
        }
    }
}
